// matrices.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <marusdefs.h>
#include <iostream>

#include <limits>

template<typename _Ty>
struct helper
{

    static inline bool    _is_null(const _Ty& value_) _THROW0()
    {
        static const _Ty eps(_STD numeric_limits<_Ty>::epsilon());
        return ((-eps <= value_) && (value_ <= eps));
    }

    static inline bool    _eq(_Ty left_,  _Ty right_) _THROW0()
    {
        return _is_null(left_ - right_);
    }

    static inline _Ty    _rectify(_Ty value_) _THROW0()
    {
        return _is_null(value_) ? _Ty(0) : value_;
    }
};

typedef size_t  size_type;


template<typename _Ty, const size_type order_>
struct refN;
template <typename _Ty>
struct ref0
{
    static void _Xran()
    {
        _THROW(_STD out_of_range, "subscript out of range in operator[].");
    }

    typedef typename ref0< _Ty > _Myt;

    _Ty&    _;
    __forceinline _Myt(_Ty* p)
        : _(*p)
    {}

    __forceinline _Myt(const _Myt& ref_)
        : _(ref_._)
    {}

    INLINE _Ty*    get() _THROW0()
    {
        return &_;
    }
    INLINE _Ty&    operator[](size_type at_)
    {
        if (at_)
            _Xran();
        return *get();
    }
};

template<typename _Ty, const size_type order_>
struct ref_select
{
    typedef typename
    _STD _If <
        order_ == 1
        , ref0<_Ty>
        , refN<_Ty, order_>
    >::type _Type;

};

template<typename _Ty, const size_type order_ = 0>
struct refN : protected ref_select< _Ty, order_-1 >::_Type
{

    typedef typename refN< _Ty, order_ > _Myt;
    typedef typename ref_select< _Ty, order_-1 >::_Type _Baset;
    _Ty&    _;

    __forceinline _Myt(_Ty* p)
        : _Baset(p)
        , _(p[order_-1])
    {}

    __forceinline _Myt(const _Myt& ref_)
        : _Baset(ref_)
        , _(ref_._)
    {}

    INLINE _Ty&    operator[](size_type at_)
    {
        if (at_ >= order_)
            _Xran();
        return get()[at_];
    }
};


template< class subscript_>
struct cell_address;

template <typename _Ty, typename Dims_>
class   _accessor;

template <const size_type size_>
struct valid_size
{
    static const size_type _size = size_;
};

template<> struct valid_size</*(const size_type)*/0> 
{
    static const int _zero_index_size_is_not_allowed = -1;
};


struct _cell
{
    typedef _cell  _subordinate;
    typedef _cell  _self;
    static const size_type _order = 0;
    static const size_type _size = 1;
    static const size_type _capacity = 1;

    typedef size_type   _path_t[1];
    typedef cell_address<_self> _address_t;

    static void _Xran()
    {
        _THROW(_STD out_of_range, "subscript out of range in operator[].");
    }

    static inline size_type _start(size_type) _THROW0()
    {
        return 0;
    }

    static inline bool _inc_index(const size_type& ) _THROW0()
    {
        return true;
    }

    static inline bool _dec_index(const size_type& ) _THROW0()
    {
        return true;
    }

    static inline bool _inc_path(const _path_t&) _THROW0()
    {
        return true;
    }

    static inline bool _dec_path(const _path_t&) _THROW0()
    {
        return true;
    }

    static inline size_type _offset(const _path_t&) _THROW0()
    {
        return 0;
    }

    //template<typename _Elem>
    //struct storage
    //{
    //    typedef typename _Elem  type;
    //};
};

template <const size_type size_, typename subordinate_ = _cell >
struct _subscript : public subordinate_
{
    typedef typename subordinate_  _subordinate;
    typedef typename _subscript<size_, _subordinate>  _self;

    static const size_type  _order = (_subordinate::_order + 1);
    static const size_type  _size = valid_size<size_>::_size;
    static const size_type  _capacity = (_size * _subordinate::_capacity);
    static const bool       _equilateral = (_size == _subordinate::_size);

    typedef size_type   _path_t[_order];
    typedef cell_address<_self>    _address_t;

    static inline size_type _start(size_type index_) _THROW0()
    {
        return index_*_subordinate::_capacity;
    }

    static inline bool _inc_index(size_type& index_) _THROW0()
    {
        if (++index_ < _size)
            return false;
        index_ = 0;
        return true;
    }

    static inline bool _dec_index(size_type& index_) _THROW0()
    {
        if (--index_ < _size)
            return false;
        index_ = _size-1;
        return true;
    }

    static inline bool _inc_path(_path_t& path_) _THROW0()
    {
        return (_subordinate::_inc_path((typename _subordinate::_path_t&)path_))? _inc_index(path_[_order-1]) : false;
    }

    static inline bool _dec_path(_path_t& path_) _THROW0()
    {
        return (_subordinate::_inc_path((typename _subordinate::_path_t&)path_))? _inc_index(path_[_order-1]) : false;
    }

    static inline size_type _offset(const _path_t& path_) _THROW0()
    {
        return (_subordinate::_capacity*path_[_order-1])+_subordinate::_offset((const typename _subordinate::_path_t&)path_);
    }

    //template<typename _Elem>
    //struct storage
    //{
    //    typedef typename _subordinate::storage::type  type[_size];
    //};
};

template< typename subscript_>
struct cell_address : public subscript_
{
    static const size_type _order = subscript_::_order;
    typename subscript_::_path_t _path;

    cell_address()
    {
        for(size_type i = 0; i < _order; ++i)
        {
            _path[i] = 0;
        }
    }

    cell_address(const cell_address& right_)
    {
        _path = right_._path;
    }

    inline
    cell_address& operator=(const cell_address& right_)
    {
        _path = right_._path;
        return *this;
    }

    inline
    size_type to_offset() const
    {
        size_type ret = subscript_::_offset(_path);
        return (ret < subscript_::_capacity)? ret : (_Xran(), -1);
    }

    inline
    bool    inc() _THROW0()
    {
        return subscript_::_inc_path(_path);
    }

    inline
    bool    dec() _THROW0()
    {
        return subscript_::_dec_path(_path);
    }

    //inline
    //cell_address& from_offset(size_type offset_) const _THROW0()
    //{
    //    return subscript_::_offset(_path);
    //}

    size_type&  operator[](const size_type order_)
    {
        if (order_ >= _order)
            _Xran();
        return _path[order_];
    }

    size_type  operator[](const size_type order_) const
    {
        if (order_ >= _order)
            _Xran();
        return _path[order_];
    }
};

template <const size_type size1_, const size_type size0_>
struct _2d : public _subscript<size0_>
{
    typedef _subscript<size0_> _subordinate;
    typedef _subscript<size1_, _subordinate>  _self;

    static const size_type  _order = _self::_order;
    static const size_type  _size = _self::_size;
    static const size_type  _capacity = _self::_capacity;
    static const bool       _equilateral = (_self::_size == _subordinate::_size);

    typedef typename _self::_address_t _address_t;

    static inline size_type _start(size_type index_) _THROW0()
    {
        return _self::_start(index_);
    }

    static inline bool _inc_index(size_type& index_) _THROW0()
    {
        return _self::_inc_index(index_);
    }

    static inline bool _dec_index(size_type& index_) _THROW0()
    {
        return _self::_dec_index(index_);
    }

    static inline bool _inc_path(_path_t& path_) _THROW0()
    {
        return _self::_inc_path(path_);
    }

    static inline bool _dec_path(_path_t& path_) _THROW0()
    {
        return _self::_dec_path(path_);
    }

    static inline size_type _offset(const typename _self::_path_t& path_) _THROW0()
    {
        return _self::_offset(path_);
    }

    //template<typename _Elem>
    //struct storage
    //{
    //    typedef typename _subordinate::storage::type  type[_size];
    //};
};

template <const size_type size2_, const size_type size1_, const size_type size0_>
struct _3d : public _2d<size1_, size0_>
{
    typedef _2d<size1_, size0_> _subordinate;

    typedef _subscript<size2_, _subordinate>  _self;

    static const size_type  _order = _self::_order;
    static const size_type  _size = _self::_size;
    static const size_type  _capacity = _self::_capacity;
    static const bool       _equilateral = ((_self::_size == _subordinate::_size) && _subordinate::_equilateral);

    typedef typename _self::_address_t _address_t;

    static inline size_type _start(size_type index_) _THROW0()
    {
        return _self::_start(index_);
    }

    static inline bool _inc_index(size_type& index_) _THROW0()
    {
        return _self::_inc_index(index_);
    }

    static inline bool _dec_index(size_type& index_) _THROW0()
    {
        return _self::_dec_index(index_);
    }

    static inline bool _inc_path(_path_t& path_) _THROW0()
    {
        return _self::_inc_path(path_);
    }

    static inline bool _dec_path(_path_t& path_) _THROW0()
    {
        return _self::_dec_path(path_);
    }

    static inline size_type _offset(const typename _self::_path_t& path_) _THROW0()
    {
        return _self::_offset(path_);
    }

    //template<typename _Elem>
    //struct storage
    //{
    //    typedef typename _subordinate::storage::type  type[_size];
    //};
};

template <const size_type size3_, const size_type size2_, const size_type size1_, const size_type size0_>
struct _4d : public _3d<size2_, size1_, size0_>
{
    typedef _3d<size2_, size1_, size0_> _subordinate;

    typedef _subscript<size3_, _subordinate>  _self;

    static const size_type  _order = _self::_order;
    static const size_type  _size = _self::_size;
    static const size_type  _capacity = _self::_capacity;
    static const bool       _equilateral = ((_self::_size == _subordinate::_size) && _subordinate::_equilateral);

    typedef typename _self::_address_t _address_t;

    static inline size_type _start(size_type index_) _THROW0()
    {
        return _self::_start(index_);
    }

    static inline bool _inc_index(size_type& index_) _THROW0()
    {
        return _self::_inc_index(index_);
    }

    static inline bool _dec_index(size_type& index_) _THROW0()
    {
        return _self::_dec_index(index_);
    }

    static inline bool _inc_path(_path_t& path_) _THROW0()
    {
        return _self::_inc_path(path_);
    }

    static inline bool _dec_path(_path_t& path_) _THROW0()
    {
        return _self::_dec_path(path_);
    }

    static inline size_type _offset(const typename _self::_path_t& path_) _THROW0()
    {
        return _self::_offset(path_);
    }

    template<typename _Elem>
    struct storage
    {
        typedef typename _subordinate::storage::type  type[_size];
    };
};

template <const size_type size_>
struct _vector 
{
   typedef typename _subscript<size_>  type;
};
template <const size_type size_>
struct _square 
{
    typedef typename _2d<size_, size_>  type;
};

template <const size_type size_>
struct _cube
{
    typedef typename _3d<size_, size_, size_>  type;
};

template <const size_type size_>
struct _hypercube
{
    typedef typename _4d<size_, size_, size_, size_>  type;
};

#include <vector>

template<typename _Elem, typename _Dim, const size_type order_>
struct storage_type_select {
    typedef typename storage_type_select< _Elem, typename _Dim::_subordinate, order_-1>::type type[_Dim::_size];
    //typedef typename
    //_STD _If <
    //    (order_ == 0)
    //    , _Elem
    //    , typename storage_type_select< _Elem, typename _Dim::_subordinate, order_-1>::type[_Dim::_size]
    //>::type type;
};

template<typename _Elem>
struct storage_type_select<_Elem, _cell, 0> {
    typedef typename _Elem type;
};


template <typename _Ty, typename _Dims>
struct  _array_base : public _Dims
{
    typedef typename _Ty _Elem;
    typedef typename _Dims::_address_t _address_t;
    typedef typename _array_base<_Ty, _Dims> _Myt;
    typedef typename storage_type_select<_Elem, _Dims, _Dims::_order>::type _storage_t;
    //typedef typename _Myt::_subordinate _Mysub;
    //typedef template _Mysub::storage<_Elem> _MysubStorage;
    //typedef typename _MysubStorage::type  _storage_t[_Dims::_size];
    //typedef typename _Dims::_subordinate::storage<_Elem>::type  _storage_t[_Dims::_size];

    typedef typename
    _STD _If <
        (_Dims::_subordinate::_order == 0)
        , _Elem&
        , _accessor<_Elem, typename _Dims::_subordinate>
    > ::type  Value;

    //typedef typename Value::storage<_Elem>::type  _storage_t[_Dims::_size];

    static inline size_type storage_size()  _THROW0()
    {
        return sizeof(_storage_t);
    }

    static inline size_type index(size_type at_)
    {
        return (at_ < _size)? at_ : (_Xran(), (size_type)-1);
    }

    static inline   size_type get_order() _THROW0()
    {
        return _order;
    }

    static inline   size_type get_size() _THROW0()
    {
        return _size;
    }

    static inline   size_type get_capacity() _THROW0()
    {
        return _capacity;
    }

    static inline   bool get_equilateral() _THROW0()
    {
        return _equilateral;
    }

    __declspec(property(get = get_order))
    size_type Order;

    __declspec(property(get = get_size))
    size_type Size;

    __declspec(property(get = get_capacity))
    size_type Capacity;

    __declspec(property(get = get_equilateral))
    bool Equilateral;
};

template <typename _Ty, typename _Dims>
class   _accessor : public _array_base<_Ty, _Dims>
{
public:
    typedef _accessor<_Elem, _Dims> _Myt;

private:
    _Myt();
protected:
    _storage_t& _storage;

public:

    _Myt(const _storage_t& storage_) _THROW0()
        : _storage(*const_cast<_storage_t*>(&storage_))
    {
    }
    _Myt(const _Myt& right_)  _THROW0()
        : _storage(right_._storage_)
    {
    }

    _Myt& operator=(const _Myt& right_)  _THROW0()
    {
        _storage(right_._storage_);
        return *this;
    }

    _Myt& operator=(const _storage_t& right_) _THROW0()
    {
        ::memcpy(memory(), right_, storage_size());
        return *this;
    }

    inline
    _Elem* memory() _THROW0()
    {
        return _storage;
    }

    inline
    const _Elem* memory() const _THROW0()
    {
        return _storage;
    }

    inline
    _storage_t& data() _THROW0()
    {
        return _storage;
    }

    inline
    const _storage_t& data() const _THROW0()
    {
        return _storage;
    }

    Value  operator[](const size_type at_) const
    {
        return (Value(_storage[index(at_)]));
    }

    _Elem    operator[](const _address_t& address_) const
    {
        return memory()[address_.to_offset()];
    }

    _Elem&    operator[](const _address_t& address_)
    {
        return memory()[address_.to_offset()];
    }
};

template <typename _Ty, typename _Dims>
class   _array : public _array_base<_Ty, _Dims>
{
public:
    typedef _array<_Elem, _Dims> _Myt;

protected:
    _storage_t  _storage;

    void initialize() _THROW0()
    {
        ::memset(memory(), 0, storage_size());
    }

public:
    _Myt() _THROW0()
    {
        initialize();
    }

    _Myt(const _storage_t& storage_) _THROW0()
    {
        ::memcpy(memory(), right_.memory(), storage_size());
    }

    _Myt(const _Myt& right_) _THROW0()
    {
        ::memcpy(memory(), right_.memory(), storage_size());
    }

    _Myt& operator=(const _Myt& right_) _THROW0()
    {
        ::memcpy(memory(), right_.memory(), storage_size());
        return *this;
    }

    _Myt& operator=(const _storage_t& right_) _THROW0()
    {
        ::memcpy(memory(), right_, storage_size());
        return *this;
    }

    inline
    _Elem* memory() _THROW0()
    {
        return (_Elem*)_storage;
    }

    inline
    const _Elem* memory() const _THROW0()
    {
        return (const _Elem*)_storage;
    }

    inline
    _storage_t& data() _THROW0()
    {
        return _storage;
    }

    inline
    const _storage_t& data() const _THROW0()
    {
        return _storage;
    }

    inline
    _storage_t& operator*()  _THROW0()
    {
        return _storage;
    }

    inline
    const _storage_t& operator*() const _THROW0()
    {
        return _storage;
    }

    Value    operator[](const size_type at_) const
    {
        return (Value(_storage[index(at_)]));
    }

    _Elem    operator[](const _address_t& address_) const
    {
        return ((_Elem*)_storage)[address_.to_offset()];
    }

    _Elem&    operator[](const _address_t& address_)
    {
        return ((_Elem*)_storage)[address_.to_offset()];
    }

};

struct _cell_coord
{
    size_type   _row;
    size_type   _col;
};

struct  round_plan_base
{
    typedef _cell_coord* cursor;

    static 
    void _reset(cursor cursor_, size_type size_ )
    {
        for(size_type i = 0; i < size_; ++i)
        {
            cursor_[i]._row = cursor_[i]._col = i;
        }
    }

    static 
    void _reset_rows(cursor cursor_, size_type size_)
    {
        for(size_type i = 0; i < size_; ++i)
        {
            cursor_[i]._row = i;
        }
    }

    static 
    void _reset_cols(cursor cursor_, size_type size_)
    {
        for(size_type i = 0; i < size_; ++i)
        {
            cursor_[i]._col = i;
        }
    }

    static
    INLINE void _swap_rows(cursor cursor_, const size_type first_, const size_type second_)
    {
        if (first_ != second_)
            _STD swap(cursor_[first_]._row, cursor_[second_]._row);
    }

    static
    INLINE void _swap_cols(cursor cursor_, const size_type first_, const size_type second_)
    {
        if (first_ != second_)
            _STD swap(cursor_[first_]._col, cursor_[second_]._col);
    }
};

template <const size_type size_>
struct round_plan : public round_plan_base
{
    typedef _cell_coord _cursor_t[size_];
    _cursor_t   _cursor;

    round_plan()
    {
        _reset(_cursor);
    }

    round_plan(const cursor right_)
    {
        memcpy(_cursor, right_, sizeof(_cell_coord)*size_);
    }

    INLINE void reset()
    {
        _reset(_cursor);
    }

    INLINE void reset_rows()
    {
        _reset_rows(_cursor);
    }

    INLINE void reset_cols()
    {
        _reset_cols(_cursor);
    }

    INLINE void swap_rows(const size_type first_, const size_type second_)
    {
        _swap_rows(_cursor, first_, last_);
    }

    INLINE void swap_cols(const size_type first_, const size_type second_)
    {
        _swap_cols(_cursor, first_, last_);
    }
};

template <typename _Ty>
class matalg
{
public:
    typedef typename _Ty _Elem;
    typedef _Elem _1x1data[1];
    typedef typename _array_base<_Elem, typename _square<2>::type>::_storage_t   _2x2data;
    typedef typename _array_base<_Elem, typename _square<3>::type>::_storage_t   _3x3data;
    typedef typename _array_base<_Elem, typename _square<4>::type>::_storage_t   _4x4data;

    class det_helper
    {
    protected:
        _STD vector<_cell_coord>  _coords;
    public:
        _cell_coord*    data()
        {
            return &_coords[0];
        }

        det_helper(size_type reserve_)
            : _coords(reserve_)
        {
            round_plan_base::_reset(data(), reserve_);
        }

    };

protected:
    const _Elem*    _data;
    const
    size_type       _stride;

public:
    matalg( const _Elem* data_, size_type stride_)
        : _data(data_)
        , _stride(stride_)
    {
    }

    ~matalg()
    {
    }

};

template <typename _Ty>
class   _Det : public matalg<_Ty>
{
public:
    typedef matalg<_Ty> _Baset;

    __declspec(noalias)
    static inline _Elem  _1x1(const _1x1data& data_) _THROW0()
    {
        return helper<_Elem>::_rectify(data_[0]);
    }

    __declspec(noalias)
    static inline _Elem  _2x2(const _2x2data& data_) _THROW0()
    {
        return helper<_Elem>::_rectify(data_[0][0] * data_[1][1] - data_[0][1] * data_[1][0]);
    }

    __declspec(noalias)
    static inline _Elem  _3x3(const _3x3data& data_) _THROW0()
    {

        return helper<_Elem>::_rectify(
               (data_[0][0] * (data_[1][1] * data_[2][2] - data_[1][2] * data_[2][1]))
             - (data_[0][1] * (data_[1][0] * data_[2][2] - data_[1][2] * data_[2][0]))
             + (data_[0][2] * (data_[1][0] * data_[2][1] - data_[1][1] * data_[2][0]))
             );
    }

    __declspec(noalias)
    static _Elem _4x4(const _4x4data& data_) _THROW0()
    {
        const _Elem* m = (const _Elem*)data_;

        _Elem ret = (m[ 0] * m[ 5] - m[ 1] * m[ 4]) // A0
                *   (m[10] * m[15] - m[11] * m[14]) // B5 
                ;

        ret -=      (m[ 0] * m[ 6] - m[ 2] * m[ 4]) // A1
                *   (m[ 9] * m[15] - m[11] * m[13]) // B4
                ;

        ret +=      (m[ 0] * m[ 7] - m[ 3] * m[ 4]) // A2
                *   (m[ 9] * m[14] - m[10] * m[13]) // B3
                ;

        ret +=      (m[ 1] * m[ 6] - m[ 2] * m[ 5]) // A3
                *   (m[ 8] * m[15] - m[11] * m[12]) // B2
                ;

        ret -=      (m[ 1] * m[ 7] - m[ 3] * m[ 5]) // A4
                *   (m[ 8] * m[14] - m[10] * m[12]) // B1
                ;

        ret +=      (m[ 2] * m[ 7] - m[ 3] * m[ 6]) // A5
                *   (m[ 8] * m[13] - m[ 9] * m[12]) // B0
                ;

        return helper<_Elem>::_rectify(ret);
    }

protected:
    __declspec(noalias)
    static INLINE _Elem  det2(const _Elem* data_, const size_type stride_, const round_plan_base::cursor plan_) _THROW0()
    {
        const _Elem* data0_ = data_ + plan_[0]._row * stride_;
        const _Elem* data1_ = data_ + plan_[1]._row * stride_;
        const size_type col0 = plan_[0]._col;
        const size_type col1 = plan_[1]._col;

        return helper<_Elem>::_rectify(data0_[col0] * data1_[col1] - data0_[col1] * data1_[col0]);
    }

    __declspec(noalias)
    static INLINE _Elem  det3(const _Elem* data_, const size_type stride_, const round_plan_base::cursor plan_) _THROW0()
    {
        const _Elem* base_row_ = data_ + plan_[0]._row * stride_;
        const round_plan_base::cursor local_plan = plan_+1;
        _Elem result(0);
        for(size_type i = 0; i < 3; ++i)
        { // 1
            _Elem val = (base_row_[plan_[i]._col]);
            round_plan_base::_swap_cols(plan_, 0, i);

            if (!helper<_Elem>::_is_null(val))
            {
                val *= det2(data_, stride_, local_plan);
                result +=  helper<_Elem>::_rectify((i&1)? -val : val);
            }
        }
        // restore columns
        round_plan_base::_swap_cols(plan_, 0, 2);
        round_plan_base::_swap_cols(plan_, 0, 1);

        return helper<_Elem>::_rectify(result);
    }

    __declspec(noalias)
    static _Elem  FASTCALL detN(size_type level_, const _Elem* data_, const size_type stride_, const round_plan_base::cursor plan_) _THROW0()
    {
        if (level_ == 1)
            return (data_+ plan_[0]._row*stride_)[plan_[0]._col];
#if 1 // ifndef _DEBUG

        if (level_ == 2)
            return det2(data_, stride_, plan_);
        if (level_ == 3)
            return det3(data_, stride_, plan_);
#endif
        const _Elem* base_row = data_ + plan_[0]._row * stride_;
        const round_plan_base::cursor local_plan = plan_+1;
        _Elem result(0);
        
        for(size_type i = 0; i < level_; ++i)
        { 
            _Elem val = (base_row[plan_[i]._col]);
            round_plan_base::_swap_cols(plan_, 0, i);

            if (!helper<_Elem>::_is_null(val))
            {
                val *= detN(level_-1, data_, stride_, local_plan);
                result += (i&1)? -val : val;
            }
        }
        while(--level_)
            round_plan_base::_swap_cols(plan_, 0, level_);
        return result;
    }

    __declspec(noalias)
    static void  FASTCALL minorN(size_type level_, const _Elem* data_, _Elem* out_, const size_type stride_, const round_plan_base::cursor plan_) _THROW0()
    {
        for(size_type i = 0; i < level_; ++i)
        { 
            _Elem* dest_row = out_ + plan_[i]._row * stride_;
            const round_plan_base::cursor local_plan = plan_+1;
            round_plan_base::_swap_rows(plan_, 0, i);
            size_type j;
            for(j = 0; j < level_; ++j)
            {
                round_plan_base::_swap_cols(plan_, 0, j);
                dest_row[plan_[0]._col] = detN(level_-1, data_, stride_, local_plan);
                //_STD cout << "row " << plan_[0]._row << ':';
                //for(size_t t = 0; t < level_; ++t)
                //{
                //    _STD cout << '['<< plan_[t]._col << ']';
                //}
                //_STD cout << _STD endl;
            }
            round_plan_base::_reset_cols(plan_, level_);
        }
    }

public:

    _Det(const _Elem* data_, size_type stride_)
        : _Baset(data_, stride_)
    {
    }

    ~_Det()
    {
    }


public:
    _Elem    det() const
    {
#if 1 //ifndef _DEBUG
        if (_stride == 1)
            return _1x1(*(_1x1data*)_data);

        if (_stride == 2)
            return _2x2(*(_2x2data*)_data);

        if (_stride == 3)
            return _3x3(*(_3x3data*)_data);

        if (_stride == 4)
            return _4x4(*(_4x4data*)_data);
#endif
        return detN(_stride, _data, _stride, det_helper(_stride).data());
    }

    void minor(_Elem* out_) const
    {
//#ifndef _DEBUG
//        if (_stride == 1)
//            return _1x1(*(_1x1data*)_data);
//
//        if (_stride == 2)
//            return _2x2(*(_2x2data*)_data);
//
//        if (_stride == 3)
//            return _3x3(*(_3x3data*)_data);
//
//        if (_stride == 4)
//            return _4x4(*(_4x4data*)_data);
//#endif
        minorN(_stride, _data, out_, _stride, det_helper(_stride).data());
    }
};


template <typename _Ty, const size_type size_>
class Vector : public _array<_Ty, typename _vector<size_>::type >
{
public:
    typedef _array<_Ty, typename _square<size_>::type > _Baset;
    typedef Vector<_Ty, size_ > _Myt;

    _Myt()
        : _array()
    {}

    _Myt(const _Myt& right_) _THROW0()
    {
        ::memcpy(memory(), right_.memory(), storage_size());
    }

    INLINE _Myt& operator=(const _Myt& right_) _THROW0()
    {
        ::memcpy(memory(), right_.memory(), storage_size());
        return *this;
    }

    INLINE _Myt& operator=(const _storage_t& right_) _THROW0()
    {
        ::memcpy(memory(), right_, storage_size());
        return *this;
    }

    static inline bool eq(const _Ty* left_, const _Ty* right_, size_type count_) _THROW0()
    {
        while(count_--)
        {
            if (!helper<_Ty>::_eq(*left_++, *right_++))
                return false;
        }
        return true;
    }

    INLINE const bool operator==(const _storage_t& right_) const  _THROW0()
    {
        return eq(memory(), right_, Capacity);//0 == ::memcmp(memory(), right_, storage_size());
    }

    INLINE const bool operator==(const _Myt& right_) const  _THROW0()
    {
        return eq(memory(), right_.memory(), Capacity); //0 == ::memcmp(memory(), right_.memory(), storage_size());
    }

};

template <typename _Ty, const size_type size_>
class square_matrix : public _array<_Ty, typename _square<size_>::type >
{
public:
    typedef _array<_Ty, typename _square<size_>::type > _Baset;
    typedef square_matrix<_Ty, size_ > _Myt;

    _Myt()
        : _array()
    {}

    _Myt(const _Myt& right_) _THROW0()
    {
        ::memcpy(memory(), right_.memory(), storage_size());
    }

    INLINE _Myt& operator=(const _Myt& right_) _THROW0()
    {
        ::memcpy(memory(), right_.memory(), storage_size());
        return *this;
    }

    INLINE _Myt& operator=(const _storage_t& right_) _THROW0()
    {
        ::memcpy(memory(), right_, storage_size());
        return *this;
    }

    static inline bool eq(const _Ty* left_, const _Ty* right_, size_type count_) _THROW0()
    {
        while(count_--)
        {
            if (!helper<_Ty>::_eq(*left_++, *right_++))
                return false;
        }
        return true;
    }

    INLINE const bool operator==(const _storage_t& right_) const  _THROW0()
    {
        return eq(memory(), right_, Capacity);//0 == ::memcmp(memory(), right_, storage_size());
    }

    INLINE const bool operator==(const _Myt& right_) const  _THROW0()
    {
        return eq(memory(), right_.memory(), Capacity); //0 == ::memcmp(memory(), right_.memory(), storage_size());
    }


    static const _Myt&  get_identity()
    {
        static 
        class _identity : public _Myt
        {
        public:
            _identity()
                : _Myt()
            {
                for(size_type i = 0; i < Size; ++i)
                {
                    _storage[i][i] = (_Ty)1;
                }
            }
        } idty;
        return idty;
    }

    __declspec(property(get = get_identity))
    const _Myt& identity;

    // returns determinant of this matrix
    _Ty det() const
    {
        return _Det<_Ty>(memory(), Size).det();
    }

    // transposes this matrix, returns reference to this matrix
    _Myt& transpose() _THROW0()
    {
        for(size_t i = 0; i < Size; ++i)
        {
            for(size_t j = i+1 ; j < Size; ++j)
            {
                _STD swap(_storage[i][j], _storage[j][i]);
            }
        }
        return *this;
    }

    _Myt& cofactorize() _THROW0()
    {
        _Elem* p = memory();
        for(size_t i = 0; i < Size; ++i)
        {
            for(size_t j = 0 ; j < Size; ++j)
            {
                *p *= (((i+j)&1)? _Elem(-1) : _Elem(1));
                ++p;
            }
        }
        return *this;
    }

    // converts this matrix to matrix of minors of right_ matrix
    _Myt& minorize(const _Myt& right_) _THROW0()
    {
        _Det<_Ty>(right_.memory(), Size).minor(memory());
        return *this;
    }

    // converts this matrix to cofactor of right_ matrix
    _Myt& cofactor(const _Myt& right_) _THROW0()
    {
        return minorize(right_).cofactorize();
    }

    // sum up the content of row
    static INLINE _Elem sum(const Value& row_) _THROW0()
    {
        const _Elem* p = row_.memory();
        _Elem ret(_Elem(0));
        for(size_type i = 0; i < row_.Size; ++i)
        {
            ret += *p++;
        }
        return ret;
    }

    // sum up the content of row
    static INLINE _Elem product_sum(const Value& left_, const Value& right_) _THROW0()
    {
        const _Elem* lp = left_.memory();
        const _Elem* rp = right_.memory();
        _Elem ret(_Elem(0));
        for(size_type i = 0; i < left_.Size; ++i)
        {
            ret += helper<_Elem>::_rectify(helper<_Elem>::_rectify(*lp++) * helper<_Elem>::_rectify(*rp++));
        }
        return ret;
    }

    // converts this matrix to its adjugate
    // when determinant is non-zero
    // otherwise does nothing
    // returns determinant of matrix
    _Elem adj() _THROW0()
    {
        _Myt tmp;
         tmp.cofactor(*this);
        _Elem ret(product_sum((*this)[0], tmp[0])); ///< calculates determinant of this matrix
        if (!helper<_Elem>::_is_null(ret))
        {
            *this = tmp;
            transpose();
        }
        return ret;
    }


    // multiplies each matrix element by scalar
    // always returns true
    bool mul(_Elem multiplier_)_THROW0()
    {
        multiplier_ = helper<_Elem>::_rectify(multiplier_);
        _Elem* mem(memory());
        for(size_type i = 0; i < Capacity; ++i)
        {
            *mem *= multiplier_; 
            if (helper<_Elem>::_is_null(*mem))
                *mem = _Elem(0);
            ++mem;
        }
        return true;
    }

    static _Elem mul(const Value& row_, const _Elem* col_) _THROW0()
    {
        _Elem result(0);
        for(size_t i = 0; i < row_.Size; ++i)
        {
            result += helper<_Elem>::_rectify(row_[i] * (*col_));
            col_ += row_.Size;
        }
        return result;
    }
    // multiplies each matrix element by scalar
    // always returns true
    _Myt& mul(const _Myt& left_, const _Myt& right_)_THROW0()
    {
        for(size_t i = 0; i < Size; ++i)
        {
            const Value& lr = left_[i];
            Value& tr = (*this)[i];

            for(size_t j = 0; j < Size; ++j)
            {
                tr[j] = mul(lr, right_.memory()+j);
            }
        }
        return *this;
    }

    // divides each matrix element by non-zero scalar
    // returns true when divisor_ != 0, otherwise does nothing and returns false
    INLINE bool div(_Elem divisor_) _THROW0()
    {
        return (helper<_Elem>::_is_null(divisor_)) ? false : mul(_Elem(1)/divisor_);
    }

    // converts this matrix to its inverse matrix
    // returns true when this matrix was successfully inverted (determinant is non-zero),
    // otherwise returns false
    INLINE bool inverse() _THROW0()
    {
        return div(adj());
    }
};
square_matrix<float, 2>::_storage_t _ =  
{
    // determinant  = 2
    { 1, 1},
    {-1, 1}
    /*
    
    */
};

// The determinant of this square_matrix is [17]
square_matrix<float, 3>::_storage_t _3 =  
{
    {1, 3, 2},
    {4, 1, 3},
    {2, 5, 2}
/* // matrix of minors is:
    {-13   2   18}
    { -4  -2   -1}
    {  7  -5  -11}
*/
};

// The determinant of this square_matrix is [24]
square_matrix<float, 4>::_storage_t _4 =  
{
    {3, 2, 0, 1},
    {4, 0, 1, 2},
    {3, 0, 2, 1},
    {9, 2, 3, 1}
};

float test_mat()
{
    const square_matrix<float, 3>& id3x3(square_matrix<float,3>().identity);
    float idDet = id3x3.det(); // must be [1]

    {
        square_matrix<float, 2> m1, m2, m3;
        m1 = _;
        m2 = _;
        m3.mul(m1,m2);
        m2.inverse();
        m3.mul(m1, m2);
        if (m3 == m1.identity)
        {
            m3.div(2);
        }
        else
            m3.div(1);
        
    }

    square_matrix<float, 3> m3;
    m3 = _3;
    square_matrix<float, 3>::_storage_t& ref=*m3;
    float ret3 = 0;
    float ret4 = 0;
    {
        ret3 = m3.det();

        square_matrix<float, 3> tmp(m3);
        tmp.inverse();
        square_matrix<float, 3>  res;
        res.mul(tmp, m3);
        res.mul(m3, tmp);
        bool iseq = (res == tmp.identity);
        if (iseq)
        {
            m3.div(2);
        }
        else
            m3.div(1);
    }   

    square_matrix<float, 4> m4;
    m4 = _4;
    {
        square_matrix<float, 4> tmp(m4);
        ret4 = m4.det();
        tmp.inverse();
        square_matrix<float, 4> res;
        res.mul(m4, tmp);
        bool iseq = (res == tmp.identity);
        if (iseq)
        {
            m4.div(2);
        }
        else
            m4.div(1);
    }
    return ret3 + ret4;
};

float test_row_col()
{
    square_matrix<float, 4> m4;
    m4 = _4;
    refN<float, 4*4> vr(m4[0].memory());
    vr[0] = 1;
    vr[3] = 4;
    vr[2] = 3;
    vr[1] = 2;
    refN<float, 4*4> vr2(vr);
    _STD cout
        << '[' << vr2[0] << ']'
        << '[' << vr2[1] << ']'
        << '[' << vr2[2] << ']'
        << '[' << vr2[3] << ']'
    << _STD endl;
    return m4.det();
}

int _tmain(int argc, _TCHAR* argv[])
{
    test_mat();
    float f = test_row_col();
	return (int)f;
}

